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The Ares I rocket uses roll control jets. These jets have aerodynamic implications as they 
impinge on the surface and protuberances of the vehicle. The jet interaction on the body can 
cause an amplification or a reduction of the rolling moment produced by the jet itself, either 
increasing the jet effectiveness or creating an adverse effect. A design of experiments test 
was planned and carried out using computation fluid dynamics, and a subsequent response 
surface analysis ensued on the available data to characterize the jet interaction across the 
ascent portion of the Ares I flight envelope. Four response surface schemes were compared 
including a single response surface covering the entire design space, separate sector 
responses that did not overlap, continuously overlapping surfaces, and recursive weighted 
response surfaces. These surfaces were evaluated on traditional statistical metrics as well as 
visual inspection. Validation of the recursive weighted response surface was performed using 
additionally available data at off-design point locations. 


Nomenclature 


A103 

= 

Ares A 103 OML, previous configuration 

A106 

= 

Ares A 106 OML, current configuration 

AOA 

= 

total angle of attack, missile axis, deg 

BTM 

= 

booster tumble motor 

CFD 

= 

computational fluid dynamics 

Ci 

= 

rolling moment coefficient 

deg 

= 

degrees 

DOE 

= 

design of experiments 

DOF 

= 

degrees of freedom 

FTV 

= 

flight test vehicle OML configuration 

E 

= 

estimate 

Ej 

= 

jet efficiency factor 

ft 

= 

foot 

JI 

= 

jet interaction 

k 

= 

number of regressors 

L 

= 

least squares operator 

LS 

= 

least squares 

lb 

= 

pounds force 

M 

= 

Mach number 

n 

= 

number of observations 

OML 

= 

outer mold line 

Patm 

= 

atmospheric pressure, lb/ft 2 

phi 

= 

roll angle, deg 

PRESS 

= 

prediction error sum of squares 

R 

= 

Rankine 

R 2 

= 

coefficient of multiple determination. 

Re 

= 

Reynold’s number 

RoCS 

= 

roll control system 
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response surface methodology 
recursive weighted sector regression 
sample standard deviation 
temperature, R 
variance 

variable inflation factors 

weighting variable 

weighting matrix 

weighting ratio 

regressor variable 

regressor matrix 

system response 

estimated system response 

total angle of attack, missile axis, deg 

least squares regressor coefficient estimate 

regressor coefficient 

change in rolling moment due to jet interaction 

change in pressure coefficient 

residual 

summation 

roll orientation, deg 

gas constant 

air density, lb/ft 2 3 

population standard deviation 

influence function 


Subscripts 

i, j,m = indexing variables 

JI = jet interaction 

RoCS = roll control system activated 


Superscripts 

T = matrix transpose 

-1 = inverse matrix 

A = estimated 
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I. Introduction 

D ESIGN of experiments (DOE) and response 
surface methodology (RSM) have been well 
understood and practiced in fields such as 
industrial engineering, automotive engineering, and 
agriculture for many decades. For the past 25 
years it has made a slow but needed infusion into 
aerospace research due in large part to the works of 
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Figure 1: Ares I major components. 


Unfortunately, it has been used very seldom in 
production cycles for math critical models of air 
vehicles where many still view the use of DOE in 
aerospace as a research project. The tests needed 
for math critical models of air vehicles consists of 
wind tunnel experiments, computational fluid 
dynamics (CFD) evaluations, and flight tests. 

These are areas of great potential for an 
improvement in information efficiency and 
accuracy due to the large amount of testing 
required. Air vehicle aerodynamics are notorious 
for being highly nonlinear, especially for body to 
body interactions, induced airflow to body 
interactions, protuberance interactions, and 
extreme flight envelope conditions such as high 
angle of attack or sideslip. The application of DOE 
and RSM to capture highly nonlinear aerodynamics 
will improve the efficiency and accuracy of these 
complex tests. 

For this research, the problem is to characterize the roll control jet interaction on the body of the Ares I rocket as 
accurately and efficiently as possible, covering the flight envelope of subsonic through supersonic speeds, and all 
possible angles of attack and roll orientations. This is a difficult problem in general to solve because it is 
experimentally very challenging to model the Reynold’s number (Re) and jet expansion physics. Therefore the 
interaction must be characterized by CFD and these results must be relied upon without validation until actual flight 
data becomes available. The CFD must use very fine grid models, 50 million grid points in this study, becoming 
computationally expensive. The Ares I flight envelope is also very expansive, requiring a large but locally accurate 
design space and model to cover all necessary conditions. Figure 1 shows a breakdown of the major vehicle 
components. 

Filling the design space can be done in several ways, one of which is to choose cardinal design Mach, angle of 
attack, and roll orientations and take a data point at each design locale. At this point there would be no modeling as 
every intersection along these cardinal conditions would contain a data point. A problem associated with this 
approach is the commonly used assumption that between the cardinal conditions the response is completely linear 
and there is no measure of the variation from linearity. If a validation point is used and the linear assumption is not 
valid, a large number of points would need to be added to capture the variation and keep the database ’’square”. On 
the other hand, when performed as response surfaces based on an orthogonal and balanced array, all linear, 
nonlinear, and factor interactions can be captured with fewer points than including all cardinal intersections. During 
testing, it is useful to evaluate the data in steps so that the test can be redirected or modified if necessary, possibly 
reducing or increasing the number of points needed, but in the fast pace of the Ares I design cycle CFD evaluation, 
this was not a possibility. In this study, the designers assumed, based on previous data, that second order designs in 
sectored design spaces would capture all physics needed for model adequacy of the jet interaction. 

Jet interaction research has been performed on jets acting normal to cylindrical bodies and flat plates 10 " 17 . This 
type of research most commonly investigates normal force and pitching moment effectiveness for abrupt control 
changes, while the current study characterizes tangential jet interactions. Jet efficiency parameters are commonly 
used as a metric to understand the effect of the jet on the body, but in this case the change in rolling moment was 
used. 
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II. Vehicle and System Description 

The Ares I rocket is the crewed launch vehicle being developed by NASA as a Space Shuttle replacement under 
the Constellation Program. It is an in-line, two stage rocket. The first stage is a five-segment reusable solid rocket 
booster that propels the vehicle to a speed of Mach 5.7 and an altitude of about 35.5 miles where it then separates 
and reenters the atmosphere. The upper stage then powers the Orion Crew Exploration Vehicle using a J-2X engine 

fueled by liquid oxygen and liquid 
hydrogen to an altitude of about 80 
miles where the Orion module 
separates and climbs to 185 miles for 
a circular orbit above Earth 18 . 

During ascent, the vehicle must 
remain within 20° of the desired 
orientation, or maintain less than 20° 
roll error. To control the vehicle roll 
within the allowed tolerance, roll 
control systems (RoCS) are used. 
The use of jet control systems is 
common on missiles and space 
vehicles. These systems have 
considerable advantages over 
conventional surface control. They 
have shorter response times and are 
effective at low dynamic pressures. 
Conventional surfaces may lack 
control authority, though the flow 
structure is very complex, and the 
control efficiency is difficult to 
predict due to the interaction of the 
jet with the external cross flow 10-17 . 

The Ares I RoCS are located near 
the rear of the interstage, forward of 
the forward frustrum as shown in 
Fig. 2. Past jet interaction studies 
have analyzed jets exiting normal to 
a body. However, the Ares I control 
jets exit tangential to the body, 
thereby attaching and wrapping 
around the body and impinging on 
protuberances. Figure 3 shows the jet 
impingement on the body and 
protuberances. Areas of warm colors 
are increased pressure on the body 
due to the jet interaction with the 
freestream flow while areas of cool 
colors represent pressure reductions. 
The jet interaction effect on the Ares 
I is very important to quantify. If the 
impact on the body causes severe 
reduction or reversal of control 
authority, a redesign may be 
necessary. 




Figure 3: CFD flow contours. 


III. Experiment Methods 

This experiment followed the guidelines for designing experiments given by Montgomery 19 and by the AIAA 
DOE focus group. These guidelines are combined below and are addressed for this research. 
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Montgomery also presents some points to keep in mind when using statistical techniques in experimentation: 

1. Use your non-statistical knowledge of the problem. 

2. Keep the design and analysis as simple as possible. 

3. Recognize the difference between practical and statistical significance. 

4. Experiments are usually iterative. 

A. Recognition and Satement of the Problem and Objectives of the Experiment. 

The RoCS jet interactions with the freestream air on the Ares I body must be properly characterized to be used in 
a simulator so that the guidance, navigation, and control engineers can evaluate the control authority of the RoCS 
system after taking the interactions into account. The characterization should be performed with as much accuracy 
and efficiency as possible, covering the flight envelope of subsonic through supersonic speeds and all possible 
angles of attack and roll orientations designated by both nominal and off-nominal trajectories during ascent. 

The nature of this experiment is characterization, not exploration or optimization. The experiment is performed 
using CFD and is known to be non-iterative because of limited CFD server time slot availability. All runs requested 
must be submitted at one time. For this reason the experiment could not be guided. 

The CFD results are deterministic in nature. A deterministic model requires that the DOE experimenter evaluate 
the necessity of adhering to the three basic principles of DOE and statistical methods. These are replication, 
randomization, and blocking. Replication is a repetition of the basic experiment and may also be referred to as 
reproducibility. Replication allows an estimate of experimental error and if the sample mean is used to estimate the 
effect of a factor, a more precise estimate of that effect can be gained. Randomization is used to force observations 
and errors to be independently distributed random variables, meaning that each point taken is not affected in any 
way by any previous data point. This assumption is required when using statistical methods. Randomization is also a 
good method for averaging out extraneous effects that are not controlled in a test environment but, in time, may 
affect the responses such as temperature changes. Blocking is used to improve the precision of response estimation 
between factors that are of no real use to the experimenter but may affect the response values. Blocking can reduce 
or eliminate variability transmitted by the presence of factors called nuisance factors 19 . For example, Favaregh and 
Fandman 7 reduced the variance in pitch damping of an electrically powered aircraft by blocking the different 
batteries used. If this experiment were empirically, not computationally based, these principles would have been 
followed while planning and executing the test matrix. However, because CFD produces no experimental error that 
can be addressed in this study, there is nothing to be gained in replicating points. Randomization is not necessary 
because there are no extraneous factors in CFD unless it is designed to produce these effects. Blocking could be 
performed with CFD if multiple grids and/or multiple individuals performed the analysis and post processing. This 
however, is out of the scope of this paper and the aero characterization problem itself but could be of interest in 
future studies. So, deterministic models such as CFD do not need to adhere to the three basic principles. 

The major gains in using DOE/RSM for deterministic experiments is the advantage of orthogonal design arrays 
to minimize the number of data points required for design space characterization and to better characterize the space. 
Orthogonality is an important property of DOE experiments in that it allows points to be added into a design space 
to capture effects while retaining uncorrelated factor coefficient estimates. In DOE terms this correlation is called 
aliasing, and it occurs when the effect of two or more terms cannot be separated into individual contributions. 

Any experiment involves risks. These could be making an inferential mistake from the data, such as concluding 
one design has lower drag compared to another when, in fact, the opposite was true, or characterizing a system 
output as being linear for a region when in fact it was highly nonlinear. For this experiment there is a risk of poorly 
characterizing the design space, missing highly nonlinear regions that cause complete loss of control, or control 
reversal problems in flight past the 20° roll error. 

B. Relevant Background 

Feveraging historical data and expert judgment can greatly benefit estimating the design space and response surface 
order required for a particular problem. If this information is available, one can decide upon the important factors 
without using screening experiments commonly used to determine the significant factors. The relevant background 
available for this data set was two databases developed using CFD and modeling methods during previous Ares I 
design cycles. The Ares I has had multiple design cycles. When these design changes exhibit significant changes 
that may affect the aerodynamics, an aerodynamic analysis must be performed. The two previous design 
configurations, which have been studied for RoCS interactions, are the FTV (also known as the Ares I-X Flight Test 
Vehicle) and the A103 (previous major configuration to the current one). The A103 and FTV have significantly 
different RoCS housing designs. This design change affects the jet interaction because of the altered pressure 
differential on the surface area of the two outer mold line (OMF) housings. RoCS housing pressurization contributes 
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greatly to the induced 
Mach rolling moment. The 

0.5 0.9 1.6 2 3 4 4.5 5 current configuration for 

this test is the A 106. An 
increase in booster 
deceleration motors (8 to 
10) on the skirt of the 
S first stage and location 

J (mid skirt to aft skirt) 

was the reasoning for a 
new analysis for A 106. 
The previous 

| databases are studied in 
all factors to determine 
the necessary 

experimental design 
choices to characterize 
Ij the response of interest 

adequately. Figure 4 
presents the database 
factor levels. Though the 
A 103 database covers a 
large area at many levels, 
the database itself was 
constructed on limited 
available CFD, and 
methods of symmetry 
were used about a roll 

orientation angle of 180°. The FTV database is comprised of CFD estimations for each database value but has only 
two angle of attack levels. Figure 5 is a typical sample of the datasets. The rolling moment is highly nonlinear 
across roll angle, nonlinear across the Mach region, and slightly nonlinear across angle of attack. 
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Figure 4: Historical databases. 



a) Roll variation. 



Figure 5: Historical jet interaction data. 


C. Choice of Factors and Ranges 

The factors for this experiment are Mach number, angle of attack, and roll orientation. These factors are 
consistent for most flight aerodynamic databases for the Ares I. The previous databases and information about the 
trajectory, such as where the maximum dynamic pressure occurs, was leveraged to select the levels and ranges for 
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this study. This study contains only Mach 0.5 to 1.6, where previous data showed the most interaction to occur. A 
constraint is determined by the convergence of the NASA Langley USM3D CFD 20 code, which does not converge 
for Mach numbers less than Mach 0.78. Therefore, the NASA Ames code OVERFLOW 21 must be used to capture 
the Mach 0.5 cases. 


Table 1: Factors, control variables. 


Factor 

Units 

Range 

Low 

Range 

High 

Priority 

Constraints 

Mach 

N/A 

0.5 

1.6 

1 

The CFD code USM3D is not 
convergent at Mach < 0.78, so the 
CFD code OVERFLOW must be 
used for Mach 0.5. 

CLj 

Degrees 

0 

8 

1 

None 

9 

Degrees 

0 

360 

1 

None 


D. Selection of the Response Variable 

For the development of this database, there was much discussion about the use of the jet efficiency factor, Ej, 
versus the jet interaction rolling moment coefficient, Jet interaction literature focuses on the jet efficiency 

terms in forces and moments 10-17 . The jet interaction rolling moment coefficient is calculated in Eq. (1) and is 
related to Ej , the jet efficiency for rolling moment in Eq. (2). For ease of implementation into the simulator, the team 
decided to use the jet interaction rolling moment coefficient, which will be the focus of the analysis. 

AClji — CirqCS ~ C l,Freestream (1) 

r? _ AC i jj+C iRoCS, Theoretical 

t J ~ r ( 2 ) 

L l.RocS, Theoretical 

E. Discussion of Design of Experiments and Response Surface Methodology 

DOE/RSM normally relies upon statistical analysis such as analysis of variance (ANOVA) using standard least 
squares (LS) modeling. This analysis technique has many well understood metrics to provide the experimenter with 
a clear understanding of the effectiveness of the data and of the extent to which the chosen candidate model 
captures the information. There are also other types of analysis, such as radial basis function modeling and Kriging 
interpolation, that can be used to provide effective results. 

DOE designs are powerfully efficient because they take advantage of the properties of orthogonal and balanced 
arrays to estimate the effects and interactions that are necessary or desired to capture the essence of the experiment. 
The experiments are guided by questions such as M How well do I need to know this information?” or ”What exactly 
am I trying to accomplish with this experiment, and how can I collect the most information in an accurate and 
efficient manner?” instead of ”How many data points can I get with this much money?” Planning an experiment 
solely on the last question might leave the analyst in a powerless position to improve the product if the goals of the 
experiment were not met and funds for additional testing are not available. Sequential analysis and building of the 
experiment is important to ensure all the goals are met during the experiment, not afterwards. Unfortunately for 
most, sequential testing is unavailable, as was the case here. The entire design was required up front with no 
opportunity for in- situ analysis. 

Each point where data is taken is referred to as an observation. This observation occurs, most likely, at set 
variable levels. These variables are called regressors, and they are normally the variables assumed to cause a 
significant change in the response. The response is the system output at the regressor levels. This response is 
commonly analyzed using least squares methods. 

Standard Least Squares. For standard least squares (LS), let n be the number of observations, in this case CFD 
points, and let k be the number of regressors (terms in the model). Assume a response variable y (rolling moment) 
with n > k observations. Each regressor variable will have a value Xy at the z th observation and y th level. The general 
model equation in terms of observations is 22 
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Vi = Po + PlXil + p 2 *i2 + "• +PkXkl 


(3) 


Po 3 " ^ ' Pj%ij 3 " 

7 = 1 

i = 1,2, ...,n 

The error terms £i are assumed to have a mean E(s) = 0, have a variance a 2 , and be uncorrelated random 
variables. The method of LS is to solve for regression coefficients, /? s, while minimizing the error. To accomplish 
this, the error term is solved for in terms of y and fis and this function is then minimized. 


n n 2 

- Po -'YjPjXij) 


i = 1 i = 1 

Here are the equations in matrix form, appropriate for a series of data. 
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To find the minimum point of error, the partial derivative of function L is taken with respect to /? and set equal to 
zero. 

= -2 X T y + 2 X T Xb = 0 (6) 

°P\ b 

X T Xb = X T y 

b = C X T XY 1 X T y (7) 

y — Xb ( 8 ) 

The solution of b is now the LS estimator of /? and y is the estimation of y. For this case, the team is not 
particularly interested in the values of b , as would be the case if trying to estimate static and dynamic derivatives of 
aircraft or rocket equations of motion. In this case the team is only interested in a good fitting model that comes 

reasonably close to the observations (y close to y), has good prediction, and has a minimum of poor characteristics 

such high variance inflation factors or low adjusted and predicted goodness of fit values (these will be explained 
further in this section). 

Weighted Least Squares. For normally distributed errors, standard LS produce estimators with good statistical 
properties. For non-normal distributions, those with longer or heavy tails, outliers in these tails may have a strong 
influence on the LS estimate, pulling the fit too much in the outlier direction. 

A single weighted LS estimation is a simple adjustment to standard LS but allows the experimenter to use 
additional knowledge of the data to produce a more robust and accurate model. Weighting is a process that explicitly 
defines the goodness of a point. For example, a time averaged observation may give differing error at different 
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factor levels, such as high a T versus low a T . If the high a T responses have a higher error, they should be assumed not 
to be as good as the low a T points and therefore, should be penalized. 

Each observation y t may have a measurement error associated with it, w t . This error may then be used as a 
weighting parameter. The experimenter shoud create a diagonal matrix of the weights. 


W = 


w 1± 

0 

0 

0 


0 0 0 

w 22 0 0 

o \ 0 

0 0 w a 


To obtain the LS estimator b , apply the weighting matrix to the least squares normal equations and estimate y 
accordingly with Eq. (8). 


b = ( X T WXY 1 X T Wy 


( 9 ) 


Useful Metrics. Statistical metrics help the analyst determine the effectiveness of the data and the adequacy of 
the regression models used. They can help determine the adequacy of the design itself and guide the experimenter in 
augmenting the design to guard against any undesired inadequacies, such as poor prediction. There are many metrics 
that the analyst may choose to guide the experiment and analysis. Some may never use any metrics, relying solely 
on visual inspection. In this analysis the prominent metrics used are variance inflation factors, predicted scaled 
variation, and goodness of fit. 

Variance inflation factors (VIFs) check for multicollinearity among the regressors. Multicollinearity is defined as 
near linear dependence among regressors that can have serious effects on model parameter estimates and the 
applicability of the model. If multicollinearity is present, the variance and covariance of the regressors can possibly 
become large 22 , depending on the level of multicollinearity. In the case that VIF values are equal to 1 , the design is 
completely orthogonal and no multicollinearity exists. In general, VIF’s should be less than 10 to indicate that 
multicollinearity is not a serious problem. In the case that the VIF is very high, one can take additional data aimed 
specifically at breaking the multicollinearity, though this may not be possible if discovered post test. Another 
method to reduce multicollinearity is to reduce the number of terms in the model. High multicollinearity does not 
necessarily translate into an inadequate model, but the individual regressor coefficients that have high VIFs will be 
poorly estimated. The model itself may be estimated well. It is important to check the design for multicollinearity if 
the goal of the experiment is to estimate the regression coefficient values. It is also important to visually inspect the 
model for over-fitting due to retaining regressors with high VIFs, which can lead to poor internal estimation and 
very poor extrapolation. For this analysis, low VIFs are certainly preferred, but overall good model fitting is a 
higher priority. 

The metric prediction variance and the property of rotatability provide measures of stability (reasonably stable 
distribution of the prediction variance 22 ) throughout the design space. Prediction variance gives an estimation of the 
predictive capability of a single point at a distance from the design origin. Rotatability is a property that indicates 
that the prediction variance of two locations in the design space that are equidistant from the design center should be 
equivalent. These two metrics rely solely on the design of the experiment and not at all on the actual responses. 

Goodness of fit metrics are well known and understood. These metrics are reasonable to use for comparing 
candidate models. The widely used R 2 value is a measure of the variability captured by the chosen regression model. 
However, using standard FS analysis, adding regression terms to the model artificially raises the R 2 value as the 
calculation is solely dependent on the sum of squares of the regression model and the total sum of squares. The more 
terms present in the model, the higher the regression sum of squares, giving no consideration of statistical 
significance. The adjusted R 2 value normalizes the goodness of fit by the number of terms in the model. So, if 
adding non-significant terms decreases the adjusted R 2 value, they should not be included. If the adjusted R 2 value 
becomes very low compared to the R 2 , the model is over defined, and either more points should be added or a lower 
order model should be used. The predicted R 2 value gives a good indication of the prediction of the model. The 
prediction error sum of squares (PRESS) is the determining factor in the predicted R 2 value. Fower PRESS is better. 
The PRESS estimates the model stability in the case that each observation was individually removed and a new 
model was fitted omitting that observation. The difference between the removed observation and the model gives a 
metric of prediction and individual observation leverage. If an observation has such high leverage that it 
dramatically influences the model, it should be investigated for validity. Replicates or an additional local inference 
space should be used to alleviate the high leverage, if possible. 
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F. Choice of DOE Experimental Design 

The present experimental design is a DOE full factorial (vertex points) augmented with axial points, constraint 
plane center points, and a mean overall center point, resulting in 23 cases for each sector, shown in Fig. 6. The full 
factorial estimates the main effects and two factor interactions. The constraint plane center and overall centers 
estimate the mean and add degrees of freedom to estimate quadratic terms. The axial points estimate quadratic 

terms. The levels were chosen so 
that all the points fell on cardinal 
Mach, roll, and a T values. The 
experiment design is repeated every 
60° of roll, which are referred to as 
individual sectors. This design is 
capable of capturing up to a full 
quadratic model with interactions 
as shown in the Taylor Series 
expansion of Eq. (10) in each 
sector. Sectors 1 and 2 share the 
points at 60° while sectors 2 and 3 
share points at 120°, and so on (Fig. 

7). 

Figure 8 presents the design 
space as a cube with the design 
points inside the cube. An 
important quality about this design 

is that it is bounding so that no extrapolation is necessary within the range of all factors. The total number of 
points to capture the nonlinear model using this design is 108 compared to 196 if a cardinal grid were used to 
capture similar effects. 

Y = + p M M + p a a+ P (p cp+ p M( JAa + p M(p M(p+p a(p a(p + Pm 2 M 2 + + P^cp 2 (10) 
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Figure 6: First sector factor levels. 
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Figure 8: Three dimensional 
view of the design space. 
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Figure 7: Design space with all sectors. 
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As mentioned in the discussion of metrics, this experiment should be analyzed for prediction stability and 
rotatability. Rotatability is retained well within individual sectors and across all sectors as shown in Fig. 9 and 10. A 
perfect model of rotatability is one with concentric circular contours. Concentric rectangles or near circular models 
are certainly adequate. Scaled prediction variance plots with completely misshapen contours, or those not centered 
about the origin, would not have the property of rotatability. Figure 9 shows the scaled prediction variance of the 
first sector, is applicable to each individual sector, and is based on the candidate model of Eq. (10). Figure 10 shows 
the full design space with the design points across all sectors. Although Fig. 10 shows some variation along the 

prediction variance contours, the general form 
of the contours shows reasonable stability. 
These prediction variance values are based on 
assuming a 6 th order model in roll, 2 nd order 
model for Mach, and a T . 



A Mach 

Figure 9: Scaled prediction variance of single sector at a T = 
4°. 


G. Performing the Experiment 

CFD modeling was used for the entire 
experiment. There is no CFD validation data for 
this particular application so the uncertainty is 
not known. Two CFD codes were used for the 
analysis. USM3D CFD analysis was performed 
at NASA Langley Research Center for Mach 
values of 0.78 to 4.5. Overflow CFD, performed 
at NASA AMES Research Laboratory, provided 
data for Mach 0.5. The external flow conditions 
were obtained from available trajectory 
information and relied solely on Mach, not on 
a T or roll angles. Table 2 shows the flow 
conditions. The jet exhaust is modeled as a 
perfect gas (y = 1.4). A RoCS firing-on case 
and a RoCS firing-off case were performed at 
each candidate point in the design space. 
Because of some communication mishaps 



across teams, not every case was performed at 
the exact factor levels as requested, but these 
were few and did not significantly affect the 
outcome of the modeling, in most cases. The 
next section discusses more on the effect of 
missing cases. 

It was assumed that data obtained at an a T = 
0° could be applied to all roll angles at that 
angle, i.e., data at a T = 0°; cp = 0° is equivalent 
to every data point at a T = 0°; cp = 0° to 360°. 
These points are not held constant in the 
response surface plots, though the final product 
would contain this feature. 


A Mach 

Figure 10: Scaled prediction variance of entire design space 
at a T = 4°. 
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Table 2: External flow conditions. 


Mach 

T (R) 

Re/ft 

P ATM (lb/ft 2 ) 

p(lb/ft 3 ) 

0.50 

514 

290,0000 

1716 

0.001935 

0.80 

491 

3,600,000 

1265 

0.001500 

0.90 

480 

3.700,000 

1113 

0.001350 

1.05 

462 

3.700,000 

907 

0.001144 

1.20 

442 

3.700,000 

738 

0.000973 

2.00 

371 

2,900,000 

277 

0.000437 

3.00 

386 

1,200,000 

80 

0.000122 

4.00 

414 

430,000 

25 

0.000034 


IV. Analysis Methods 


A. Analysis of the Data 

Four types of analysis were performed for evaluation. First, an analysis fitting one high order polynomial to the 
entire data set was performed. Second, quadratic models were fitted to individual sectors, and the estimated 
responses were averaged where the surfaces join. Third, a quadratic model was fit in 60° roll sectors every 15°, 
creating constantly overlapping response surfaces that can be averaged in the entire design space. Fourth, 
overlapping sectors from the third analysis were used along with recursive weighting techniques to minimize the 
maximum roll coefficient derivative when retaining CFD values. 



Mach Roll 

Figure 11: Jet interaction rolling moment vs all factors. 

Plotting the response versus all factors before fitting the model can give the analyst an idea of the main effects 
before getting into the analysis. Figure 1 1 contains three subplots that contain all the response values in each plot 
with an averaging line between groups of points. Clearly the response generally decreases as Mach increases. A 
correlation of Mach and the difference in total pressure and jet exit pressure may exist, but was not considered in 
this analysis. The response shows a decrease with a T but not as strong as the Mach effect. The response in roll shows 
no trend, indicating either no dependence or highly nonlinear behavior that is dependent on interactions with the 
other factors. The maximum value across roll is the repeated a T = 0°, cp = 0° case. 

In the following sections a response surface of the estimated model and a response surface with CFD are plotted. 
When plotting the estimated model alone, the response is of the general form y = Xb, herein referred to as Model A. 
When including the design points, the response is of the general form y = Xb + £, herein referred to as Model B. 
The distinction highlights an important discussion for deterministic experiment model analysis. Should the absolute 
values remain in the model or should the general trends of the model stand alone? A pro for adding the residual, £, 
back to the model is that at those points in the model, the estimation is exact, and since there is no replication for 
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deterministic models, there should not exist another value for that locale. However, the prediction of new points near 
that locale will not necessarily be made better by including those points. In fact, the newly predicted points may be 
worse if the residual has a high leverage on the model. The confidence intervals of the model will also be affected. If 
the model includes the residuals, then no residuals exist for a residual analysis, therefore the only means of obtaining 
error is to take validation points. While taking validation points is always a good idea to be sure the model is 
predicting well in the design space, it should not be the only means of evaluating prediction error. It is the author’s 
opinion that if the model can be made smooth with the addition of the residuals to the model, that model B should be 
used and the error from the model A should be applied to model B. 

B. Single Regression Model Analysis 

The single regression model has degrees of freedom (DOF) available for 2 nd order Mach and angle of attack 
terms and enough DOF for a 12 th order roll term. The goodness of fit is optimal, however, at the 5 th order roll term. 

The model structure for this analysis is: 

AQj/ = Po + PmM + Pa a + PcpV + PmccMoc + p M(p Mcp + (3 a( pCC(p + /?m 2 M 2 + (3 a2 a 2 + P(p 2 (p 2 + P(p 3 <P 3 + 

+ (“) 

Overall the model performs very well. Figure 12 presents the single regression model at a T = 4°. The model 
performs worst at a Mach number of 1.32 and a T = 2°. In this region, the data has high amplitude that is not 
captured well by this global model. Figure 13 includes the design points into the response at Mach 1.32 and clearly 
shows how the CFD, when added back in, largely pulls away from the surface, as well as creates large derivative 
changes where there exist two CFD points and one regression model point between them. This indicates that there 
is still something to be gained over the single regression model. The model is hierarchal, i.e., it includes all terms of 
lower order. The model goodness of fit (R 2 ) values are R 2 = 0.9488, adjusted R 2 = 0.9431, and predicted R 2 = 
0.9348. All goodness of fit values indicate that most of the variance in the data is captured by the model. The R 2 is 
reasonably high, though a higher value is always desirable. The adjusted R 2 is well in line with the R 2 , indicating 
that the model is not being over- fitted to capture the data well and that 94.31% of the variability in the data is 
captured by the model, leaving almost 6 % to fitting error. A predicted R 2 value within 2% of the R 2 suggests very 
good robustness in the model and insensitivity to missing data or highly leveraged points. The VIFs are very high 
(>10) for higher order roll terms. In order to reduce them below a value of 10, the model can be reduced to 3 rd order, 
though when reduced, the model does not capture the desired trends in roll. This is an attribute of the highly 
nonlinear nature of the aerodynamic system. 



Mach 

Figure 12: Single regression model. 


C. Multiple Sector Regression 
Analysis 

The initial plan of the 
experiment design was to capture 
up to 2 nd order behavior every 60° 
of roll. For the multiple sector 
analysis, backward elimination 
least squares (LS) regression was 
used to develop analytical models 
for each sector. Each sector is 
estimated initially by the model in 
Eq. (10) and insignificant terms 
removed. This method produced a 
combination of reduced and full 
models, depending on the 
significant variables in the sector 
spaces. Figure 14 shows the 
goodness of fit values. 

The majority of R 2 , adjusted R 2 , 
and predicted R 2 values seems very 
reasonable and indicates adequate 
modeling of the variance in the 
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Jl Rolling Moment 



Figure 13: Single regression model including design points. 
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Figure 14: Goodness of fit metrics for multiple sectors. 



Figure 15: Multiple sector response surfaces. 


data. The exception is sector 2, 
which has a significantly low 
predicted R 2 value. Figure 15 
shows the multiple response 
surfaces with the corresponding 
CFD points (not included in the 
model) at a T = 4°. Each sector 
modeled the data well, 
individually, with the exception of 
the 60° to 120° sector (sector 2). 
Unfortunately, sector 2 is missing a 
data point at Mach 0.5, a T = 8°, cp = 
120°. The sector was not bounded 
as was designed, extrapolation was 
necessary, and, therefore, an 
extremely poor prediction ensued, 
as shown by the low predicted 
R 2 values and visual inspection of 
the response surface with the 
corresponding data. 

The high sensitivity of this 
approach to missing data is 
highlighted with the red tipped 
spike at 120° roll. This poorly 
predicting response is not present 
in the single regression model 
analysis. Significant differences 
also exist in response models 
where they join. These differences 
will cause a final averaging of the 
sectors to contain undesirably high 
derivative changes across roll. The 
averaged R 2 value for the global 
model = 0.9811, adjusted R 2 = 
0.9698, and predicted R 2 = 0.9381. 
The VIFs were very low for most 
sectors, indicating a near 
orthogonal design for the chosen 
modeling as was planned, again 
with the exception of sector two, 
which had some collinearity 
because of the missing data. 

While this analysis was the 
planned choice for modeling the 
data, and the overall goodness of fit 
values show a better fit than the 
single response surface analysis, 
the consequences of this modeling 
are unacceptable. The author then 
attempted another type of modeling 
to alleviate the pitfalls of this one 
by applying an extension of the 
multiple sector analysis. 
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D. Continuous Sector Regression Analysis 

The aim of continuous sector regression analysis is to create a globally robust model, insensitive to missing data, 
and highly locally accurate in regions where the single regression model was inadequate. The models are built in 60° 
roll spaces, overlapping by 15°. The model used is the same as Eq. (10). This scheme allows each point in the design 
space to be estimated by four response surfaces. Each estimate is then averaged to create an averaged global surface. 
Advantage is taken in the circular referencing of roll; i.e., 15° = 375° and 30° = 390° etc. This property ensures all 
points are covered by four response surfaces. 

Before using the original test matrix 
in an unplanned fashion, the matrix 
must be tested for adequate properties 
in these in-between sectors. Figure 16 
represents the available data in the first 
offset sector from roll values of I5to 
75°. Figure 17 shows that the data is 
very near rotatable, a desired attribute. 
The VIF are very low for a 2 nd order 
model within this design space. It is 
feasible then to fit a 2 nd order model to 
this data in each offset sector. A 
downside to the offset sector design 
space is that the model is extrapolating 
for values of Mach 0.5 and 1.6 at roll 
angles of 15° and 75°. Fortunately, there 
are four estimations at each point, with 
only one as an extrapolated value, so 
the extrapolation effect is minimal. 

Figure 18 shows the continuous 
sector response surface. From both 
visual inspection and model quality 
metrics, the model seems to perform 
well. When retaining the CFD values in 
the response surface, it blends 
smoothly, with the exception, again, at 
Mach 1.32, a T = 2, though in this 
region it does perform better than the 
single regression model. The overall 
goodness of fit values are R 2 = 0.965, 
adjusted R 2 = 0.960, and the predicted 
R 2 = 0.964. A predicted R 2 value that is 
almost identical to the R 2 value is a 
clear sign of robustness in the 
modeling, meaning that even though 
data was not available for some cases, it 
has a negligible effect on the resulting 
model, unlike the multiple sector 
analysis. 


Figure 16: First offset sector available data points (not pure design 
points). 



A: MACH 

Figure 17: Scaled predicted variance for the first offset 
sector. 
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Mach 

a) Without design points. 



Math 


b) With design points. 


Tt J 



c) With design points at Mach 1.32. 


Figure 18: Continuous sectors. 
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E. Recursive Weighted Sector Regression Analysis 

The goal of performing recursive weighted sector regression (RWSR) analysis is to increase the accuracy and 
smoothness of local regression models. This is attempted by first obtaining a global model, then using estimations 
from the global response surface as down weighted system response points, while using the true CFD observations 
as highly weighted values. The global model could be that from the single regression model or from the continuous 
sector model because each has similarly good properties. The weighted LS Eq. (9) was used to calculate new LS 
estimators and y = Xb + £ obtained the new model estimate. A code was developed to perform this iteration on 
mini-sectors using roll increments of only 5° and performing the RWSR on 30° wide sectors (in roll). The CFD and 
response surface estimations were weighted using the variance associated with the CFD response. This creates an 
equal weighting that must then be varied using a weighting ratio, shown in Eq. (12), where w is the diagonal weight 
matrix. The reason for showing the denominator in this equation, even though it has no effect in this case, is that it 
creates a placeholder for using differing weighting denominators, such as those that can be used if wind tunnel and 
CFD results were combined. The CFD weight is varied to understand the effect of different weighting ratios ( WR ). 


w(ii) = 


var(y ) ' 


if RS 


WR 

var(y) ’ 


if CFD 


( 12 ) 


The maximum response roll derivative was tracked until a minimum was reached to ensure smoothness in roll, 
shown in Fig. 19. Examining the WR in the response surface metrics gave an insight into what weighting can do to a 
surface. The higher the weight given to the CFD, the more the surface is pulled to those points. This has a direct 
effect on the residuals of the true data to the surface. The higher the WR values, the lower the residuals. Also, as the 
WR is increased, the maximum jet interaction (JI) rolling moment derivative with respect to roll was minimized 
exponentially faster, and the value of the derivative was lower for higher WR. 



1 - 


0 1 J ! 
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Figure 19: Maximum JI rolling moment derivative with 
respect to roll when including the design points. 



Weighting Ratio (CFD/RS) 


The drawback of high WR is the 
considerably high flexibility in the model. This 
high flexibility must be restrained to reasonable 
levels to not induce too much variation in the 
model in order to reach the outer lying points. 
The influence of the WR on goodness of fit is 
shown on a log scale in Fig. 20. 

The R 2 values increase logarithmically, 
implicitly depicting the decrease in residuals. 
So the art is to select the WR to allow smooth 
transition to true CFD points without inducing 
undesired variability and retaining a low 
maximum derivative. In this case a value of 10 
was chosen that reduced the maximum 
derivative 30% (Fig. 19). 

The overall goodness of fit values are R 2 = 
0.9806, adjusted R 2 = 0.9782, and the predicted 
R 2 = 0.9804. Although the goodness of fit 
value gain seems to be very marginal over that 
of the previous methods, the availability to 
include the deterministic points appropriately 
makes this approach the most robust (it is least 
influenced by missing data) and most accurate 
(it captures the known data well while not 
inducing unreasonable variation in the model). 
This approach is also simple to program and 
implement. Figure 21 presents the final model 
with WR = 10 at Mach 1.05. 


Figure 20: Influence of weighting ratios. 
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Figure 21: RWSR model at Mach 1.05 including design points. 


Montgomery 22 presents an 

optimization technique called iteratively 
reweighted least squares that can be 
applied to this problem of including 
deterministic data into a surface smoothly 
as well as blending deterministic and 
experimental results. The application of 
this method is worthwhile for future work 
and is briefly introduced. 

While the name sounds similar to the 
recursive weighting method just discussed, 
the weighting is different. To perform 
iteratively reweighted least squares, an 
initial estimation of b is obtained, indexed 
as b 0 . The weighting is then a diagonal 
matrix with elements w i0 such that 


w, 


io 


(ip[(yi - x[b 0 )/s], if y i * x[b 0 
[ 1 ,ify i =x' i b 0 


(Ref. 22) 


(13) 


b m = (X T W m xr 1 X T W m y 


(14) 


where s is a robust estimate of the standard deviation, i is the observation, m is the iteration index, and x/j is an 
influence function that controls the individual weights of residuals, applicable within a range of tuning factors. The 
function of residuals, which is to be minimized, is called a robust criterion function. Montgomery 22 gives a handful 
of commonly used robust criterion functions with corresponding influence functions and weights. Eq. (9) is 
modified to provide iterative coefficient estimations. 


F. Validation 
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Figure 22: Residuals to model error and validation. 


No model should be relied upon without having estimated error predicting points inside the area of operation. 
Validation points were used to evaluate the predictability of the model at off-design points within the area of 
operation of this response surface. In a wind tunnel test setting, these points would be used to exercise the model and 
guide the experiment into regions if the prediction is not within tolerance. For this CFD study, validation directly 
aids the development of an uncertainty model for the response surface. The RWSR model was used to predict the 
response at 28 available validation points 1 . These points were not used in the creation of the response surface. The 
validation points were at Mach 0.9 and 1.2. Residuals of the RWSR predictions and the validation residuals are 
plotted in Fig. 22. A three standard deviation dispersion bound is plotted for the model-only residuals and for the 


^ Even with the additional 28 validation points, the total for this experiment is 60 fewer points than for the cardinal grid. Each 
case consists of a power on and power off case, so the real savings is 120 CFD runs. 
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model + validation residuals. While these bounds do not represent the actual bounds that would be chosen given this 
data set, they do show the general fitness of using bounds created without validation points compared to bounds 
created with them. Mach 0.9 validation points show deviation from the model, which is not surprising given the flow 
physics associated with transonic flow near Mach 1. An inspection of the residuals versus Mach indicates an 
uncaptured 3 rd order Mach variation. The model could be refitted including the validation points with a (3 M3 M 3 term 
to create a more constant variation across Mach. The residuals are mostly constant across a T , but there is certainly 
more variation in the second sector of roll, 60 to 120°, than in the rest of the sectors. If experiment iteration were 
possible, an augmentation to the design in this sector would have added the additional DOF to capture the large 
variations in roll. Having constant error across all factors is preferred. It means that the model is going through the 
mean of the data at all levels and that all variance from factors is captured. If the error is constant, the uncertainty 
becomes a simple value that can be applied across the entire space. 

Figure 23 shows the modeling residuals when including the available validation points in the modeling. These 
plots show a constant variance across all variables; no longer are the residuals varying with Mach and roll. There is a 
significant improvement over the model constructed without the validation cases. Had sequential testing been 
performed, these cases would have proven very valuable to guiding the test. Luckily, having the foresight to check at 
these conditions turned out to be a sound decision. 



0.6 0.8 1 1.2 1.4 1.6 0 2 4 6 8 0 30 60 90 120 150 180 210 240 270 300 330 360 
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Figure 23: Residuals with validation points included during modeling. 

Referring back to the discussion of Model A and Model B, B including the observation values and A not 
including them, the final product for this dataset would most likely include the validation points while retaining the 
error associated with predicting the validation points. This would include all relevant information and retain a good 
estimate of prediction error within the area of operation. This method allows for constant and easy model and error 
updating when new data becomes available. 

G. Lessons Learned and Future Work 

All values at a T = 0° were applied across roll, i.e., the value at (p = 0° is the constant value across all values of (p 
for that Mach number. This seems to gain efficiency in the data matrix because only one CFD run is taken for all roll 
values at zero a T . But, this creates very large derivative changes between a T = 0° and a T = 2°, where the response is 
very nonlinear. This is problematic because the surface has to ’’catch up” to the nonlinearities or smooth to a 
constant value at a T = 0°. In future studies, a T = 0° values will be taken in addition to values at incidence, as a 
separate design. Because so few are required, it will not impact the efficiency of the design. The next design should 
contain an incidence range of 2 to 8°. Much diligence should be taken to ensure that all points in the design space 
are taken at the exact values as requested, in order to guard against multicollinearity and protect the prediction 
capabilities of the response surfaces. 

Jet efficiency should be evaluated, in addition to the rolling moment variable in this study, in the same fashion. 
The efficiency factor could be much less complex and require less complex modeling. If this is the case, a sequential 
test schedule, using efficiency factor as the main response, could result in a dramatically reduced test matrix. 

The use of non-parametric modeling algorithms should be evaluated and compared to the results of this study. 
Methods such as Kriging interpolation and radial basis functions were not available for this work but could have 
proven more successful than the LS methods studied here. 
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V. Conclusion 

The application of DOE/RSM to highly nonlinear aerodynamics characterization problems can be done well 
using some of the techniques presented in this work. All models except for the multiple sector, characterized the 
data adequately. The downfall for the multiple sector analysis was the missing bounding data points in sector 2, 
highlighting the sensitivity of the approach to missing data. This approach could still be employed in the future if all 
requested observations are available, as most individual sectors performed very well. The goodness of fit values of 
each approach are summarized in Table 3. It is apparent from this comparison that the multiple sector approach 
performed almost as well as the recursive weighted approach except for the prediction R 2 values. The single 
regression model is the easiest to perform and evaluate. This model can also be optimized more readily as there is 
only one function to define it. In cases where model accuracy of 93-95% is acceptable, this would certainly be the 
desired approach. The continuous sector approach performed adequately, but the addition of complexity did not gain 
much benefit over the single model. However, combining the continuous sector method with recursive weighting 
and response surface observations provided the most accurate and usable models, as the inclusion of true values did 
not induce large derivative changes that are undesired. This approach has promise for use in wind tunnel testing as 
well; where individual point goodness can be used to weigh replicates. The model can be updated in-situ as each 
new point is taken until a tolerance in prediction or some other factor has been obtained. A similar type of in-situ 
dynamic modeling has been applied effectively and accurately to flight test data by Morelli and Smith 23 . An 
extension to ground testing should be natural. The application of robust criterion and influence functions will also 
expand the application of DOE/RSM to highly nonlinear systems. 


Table 3: Goodness of fits. 


Goodness of Fit 

Single 

Multiple 

Continuous 

Recursive Weighted 

R 2 

0.949 

0.981 

0.965 

0.981 

Adjusted R 2 

0.943 

0.970 

0.960 

0.978 

Prediction R 2 

0.935 

0.938 

0.964 

0.980 


The application of these types of experiment designs and modeling schemes can be used in deterministic 
systems as well as those with systematic and random errors. The gain in efficiency of orthogonal arrays can be 
effective in reducing design cycle time and resource needs while gaining more information per individual 
observation than can be accomplished in traditional test methods. 
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